clear all
set more off

global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"


cap log close
log using $figures\N_file,replace t

**************************************************************************************************************
cap program drop pdm1m0
program pdm1m0
	args x1 x2 r
	g pDi=binormal(`x1',`x2',`r')/normal(`x2')						/*This is Pr(D=1|A=1)*/
	g m1i=normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',`x1',`r') + ///
		`r'*normalden(`x2')*(normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',`x1',`r')
	gen m0i=-normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',-`x1',-`r') + ///
		`r'*normalden(`x2')*(1-normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',-`x1',-`r')	
	scalar sigma=sqrt(1+s_csi^2^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)

	qui su part1 if F==1 & A==1
	scalar p1_1=r(mean)
	qui su part1 if F==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	qui su part2 if F==1 & A==1
	scalar p2_1=r(mean)
	qui su part2 if F==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	qui su part3 if F==1 & A==1
	scalar part3m=r(mean)
	
	qui su pDi if F==1 & A==1	
	scalar pD_1=r(mean)
	qui su pDi if F==0 & A==1	
	scalar pD_0=r(mean)
	qui su m1i if F==1 & A==1	
	scalar m1_1=r(mean)
	qui su m1i if F==0 & A==1	
	scalar m1_0=r(mean)
	qui su m0i if F==1 & A==1	
	scalar m0_1=r(mean)
	qui su m0i if F==0 & A==1	
	scalar m0_0=r(mean)
	
	scalar drop sigma p1_1 p1_0 p2_1 p2_0
	drop pDi m1i m0i part1 part2 part3
end
**************************************************************************************************************

*********************************************************************************************************************************
/*CORRECTION FACTORS*/
cap program drop get_correction
program get_correction
	scalar cfl   =bF_L-((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2))) 
	scalar cfa   =bF_A-((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))
	scalar cfr   =bF_R-delta
	scalar cfcla =rho_LA-(sqrt(1+s_omega^2+s_psi^2)/sqrt(1+s_omega^2+s_psi^2+s_v^2))
	scalar cfclr =rho_LR-(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2)))
	scalar cfcra =rho_RA-(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2+s_v^2)))
end
***********************************************************************************************************************************

cap program drop get_baseline_pars
program get_baseline_pars
	preserve
		u $data/R_estimates_DWMD,clear
		destring DWMD,gen(beta)
		replace namev="alpha" if namev=="alpha1"
		ren namev namepar
		*********************************************************************************************************************************
		*BASELINE PARAMETERS
		*********************************************************************************************************************************
		su beta if namepar=="pi"
		scalar mu		=r(mean) 
		su beta if namepar=="gamma"
		scalar gamma	=r(mean)
		su beta if namepar=="tau"
		scalar tau		=r(mean)
		su beta if namepar=="sigma_omega"
		scalar s_omega	=r(mean)
		su beta if namepar=="sigma_psi"
		scalar s_psi	=r(mean)
		su beta if namepar=="sigma_v"
		scalar s_v		=r(mean)
		su beta if namepar=="sigma_csi"
		scalar s_csi	=r(mean)
		su beta if namepar=="alpha"
		scalar alpha1	=r(mean)
		scalar sigma	=sqrt(1+s_csi^2)			/*This is [1+var(noise_SSA)] */	
		***********************************************************************************************************************************
		*********************************************************************************************************************************
	restore
end


cd $data
***********************************************************************************************************************************************************************************
u $data\after_ML, clear

scalar rho_LA=rho_A_L
scalar rho_RA=rho_A_R
scalar rho_LR=rho_L_R
scalar betaL=bF_L
scalar betaA=bF_A
scalar betaR=bF_R

ren female F
gen Aw=1-R
ren logbs Y

ren a1 Lhat					/*This includes the effect of gender*/
ren a2 Ahat					/*This includes the effect of gender*/		
ren a3 Rhat					/*This includes the effect of gender*/

gen xbL=Lhat-bF_L*F			/*subtract the effect of gender because it is going to add it "structurally" later*/
gen xbA=Ahat-bF_A*F			/*subtract the effect of gender because it is going to add it "structurally" later*/
gen xbR=Rhat-bF_R*F			/*subtract the effect of gender because it is going to add it "structurally" later*/	


get_baseline_pars

gen Dhat=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat Ahat rho_DA

scalar delta=part1m+part2m+part3m*alpha1 	

get_correction

gen cfid=1
replace cfid=sum(cfid)
				
matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)
putmata cfid U1=(Lhat Ahat Rhat) if !missing(Rhat), replace

mata
	R1=vech(st_matrix("W1"))'
	Z=mvnormal(U1,R1)
	st_matrix("Z",Z)
end

getmata Z, id(cfid) replace
g Type1_baseline=Z/binormal(Ahat,Lhat,rho_LA)
su Type1* if F==1 & A==1
su Type1* if F==0 & A==1
save trial_figure1,replace


******pi*******
get_baseline_pars

matrix define mu_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar mu_x=`j'

	g Lhat_pi=xbL+(((mu_x-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_pi=xbA+(((mu_x-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_pi=/*xbLV*/xbL*sqrt(1+s_omega^2+s_psi^2)+mu_x*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_pi Ahat_pi rho_DA
	scalar delta=part1m+part2m+part3m*alpha1
	g Rhat_pi=xbR+(delta+cfr)*F
	matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)			
	putmata cfid U1=(Lhat_pi Ahat_pi Rhat_pi) if !missing(Rhat_pi), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type1_pi=Z/binormal(Ahat_pi,Lhat_pi,rho_LA)

	qui su Type1_pi if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type1_pi if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix mu_results=mu_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat mu_results
keep if mu_results5!=.
keep mu_result*
drop if mu_results5==0
ren mu_results2 	pi_Type1_F
ren mu_results3		pi_Type1_M
ren mu_results4		pi_m1_F
ren mu_results5		pi_m1_M
ren mu_results6		pi_m0_F
ren mu_results7		pi_m0_M
ren mu_results8		pi_pD_F
ren mu_results9		pi_pD_M
ren mu_results1 	value
sort value
save mu_results,replace


******gamma*******
get_baseline_pars

matrix define gamma_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar gamma_x=`j'

	g Lhat_gamma=xbL+(((mu-gamma_x)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_gamma=xbA+(((mu-gamma_x-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_gamma=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_gamma Ahat_gamma rho_DA
	scalar delta=part1m+part2m+part3m*alpha1
	g Rhat_gamma=xbR+(delta+cfr)*F
	matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)			
	putmata cfid U1=(Lhat_gamma Ahat_gamma Rhat_gamma) if !missing(Rhat_gamma), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type1_gamma=Z/binormal(Ahat_gamma,Lhat_gamma,rho_LA)

	qui su Type1_gamma if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type1_gamma if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix gamma_results=gamma_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat 	gamma_results
keep if gamma_results5!=.
keep 	gamma_result*
drop if gamma_results5==0
ren gamma_results2 		gamma_Type1_F
ren gamma_results3		gamma_Type1_M
ren gamma_results4		gamma_m1_F
ren gamma_results5		gamma_m1_M
ren gamma_results6		gamma_m0_F
ren gamma_results7		gamma_m0_M
ren gamma_results8		gamma_pD_F
ren gamma_results9		gamma_pD_M
ren gamma_results1 	value
sort value
save gamma_results,replace


******tau*******
get_baseline_pars

matrix define tau_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar tau_x=`j'
	g Lhat_tau=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_tau=xbA+(((mu-gamma-tau_x)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_tau=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_tau Ahat_tau rho_DA
	scalar delta=part1m+part2m+part3m*alpha1
	g Rhat_tau=xbR+(delta+cfr)*F
	matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)			
	putmata cfid U1=(Lhat_tau Ahat_tau Rhat_tau) if !missing(Rhat_tau), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type1_tau=Z/binormal(Ahat_tau,Lhat_tau,rho_LA)

	qui su Type1_tau if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type1_tau if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix tau_results=tau_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat 	tau_results
keep if tau_results5!=.
keep 	tau_result*
drop if tau_results5==0
ren tau_results2 		tau_Type1_F
ren tau_results3		tau_Type1_M
ren tau_results4		tau_m1_F
ren tau_results5		tau_m1_M
ren tau_results6		tau_m0_F
ren tau_results7		tau_m0_M
ren tau_results8		tau_pD_F
ren tau_results9		tau_pD_M
ren tau_results1 	value
sort value
save tau_results,replace


******alpha*******
get_baseline_pars

matrix define alpha_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar alpha_x=`j'
	g Lhat_alpha=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_alpha=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_alpha=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_alpha Ahat_alpha rho_DA
	scalar delta=part1m+part2m+part3m*alpha_x
	g Rhat_alpha=xbR+(delta+cfr)*F
	matrix W1=(1, rho_LA, rho_LR \ rho_LA,1,rho_RA \ rho_LR, rho_RA, 1)			
	putmata cfid U1=(Lhat_alpha Ahat_alpha Rhat_alpha) if !missing(Rhat_alpha), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type1_alpha=Z/binormal(Ahat_alpha,Lhat_alpha,rho_LA)

	qui su Type1_alpha if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type1_alpha if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix alpha_results=alpha_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat 	alpha_results
keep if alpha_results5!=.
keep 	alpha_result*
drop if alpha_results5==0
ren alpha_results2 		alpha_Type1_F
ren alpha_results3		alpha_Type1_M
ren alpha_results4		alpha_m1_F
ren alpha_results5		alpha_m1_M
ren alpha_results6		alpha_m0_F
ren alpha_results7		alpha_m0_M
ren alpha_results8		alpha_pD_F
ren alpha_results9		alpha_pD_M
ren alpha_results1 	value
sort value
save alpha_results,replace




u mu_results,clear
merge value using gamma_results
drop _merge
sort value
merge value using tau_results
drop _merge
sort value
merge value using alpha_results
drop _merge
sort value


foreach x in Type1_F Type1_M m1_F m1_M m0_F m0_M pD_F pD_M {
	lab var pi_`x' "{&pi}"
	lab var gamma_`x' "{&gamma}"
	lab var tau_`x' "{&tau}"
	lab var alpha_`x' "{&alpha}"
	}
lab var value "Parameter value"

scatter *_Type1_F value,c(l l l l) s(i o dh +) clp(dash solid solid solid) legend(col(4) pos(6)) graphr(c(white)) ytitle(Pr(Type I error)) ///
	xlabel(-1(.2)1) 
graph export $figures\figureOA2.pdf,replace
	



**********************************************************************************
**********************************************************************************
**********************************************************************************
**********************************************************************************
**************************************************************************************************************
**************************************************************************************************************
cap program drop pdm1m0
program pdm1m0
	args x1 x2 r
	g pDi=binormal(`x1',`x2',`r')/normal(`x2')						/*This is Pr(D=1|A=1)*/
	g m1i=normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',`x1',`r') + ///
		`r'*normalden(`x2')*(normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',`x1',`r')
	gen m0i=-normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',-`x1',-`r') + ///
		`r'*normalden(`x2')*(1-normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',-`x1',-`r')	

	scalar sigma=sqrt(1+s_csi^2^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)

	qui su part1 if F==1 & A==1
	scalar p1_1=r(mean)
	qui su part1 if F==0 & A==1
	scalar p1_0=r(mean)
	qui su part2 if F==1 & A==1
	scalar p2_1=r(mean)
	qui su part2 if F==0 & A==1
	scalar p2_0=r(mean)
	qui su part3 if F==1 & A==1
	scalar p3_1=r(mean)
	qui su part3 if F==0 & A==1
	scalar p3_0=r(mean)
	drop part1 part2 part3 pDi m1i m0i
	scalar drop sigma
end
**************************************************************************************************************
**************************************************************************************************************
matrix define c1c2_results=(0,0,0)

forvalues j=0.01(0.01)8	{
	u trial_figure1,clear

	scalar c1_c2=`j'

	g Lhat_pi=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_pi=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_pi=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

	pdm1m0 Dhat_pi Ahat_pi rho_DA

	scalar delta_1=p1_1+p2_1+p3_1*(c1_c2)
	scalar delta_0=p1_0+p2_0+p3_0*(c1_c2)
	
	qui su xbR if F==1 & A==1
	scalar xbR_1=r(mean)
	qui su xbR if F==0 & A==1
	scalar xbR_0=r(mean)
	
	g award_1=1-normal(xbR_1+(delta_1+cfr))
	g award_0=1-normal(xbR_0+(delta_0+cfr))

	qui su award_1 if F==1 & A==1			
	scalar awf=r(mean)
	qui su award_0 if F==0 & A==1
	scalar awm=r(mean)
	*********************************************************************************************************************************
	matrix c1c2_results=c1c2_results\(`j',awf,awm)
	*********************************************************************************************************************************
	}

g pDi=binormal(Dhat_pi,Ahat_pi,rho_DA)/normal(Ahat_pi)		/*it's independent of alpha, so anyone is ok*/
su pDi if F==1 & A==1
scalar pD_F=r(mean)
su pDi if F==0 & A==1
scalar pD_M=r(mean)
global pdm = pD_M
global pdf = pD_F
global arm = 0.46	/*From paper*/
global arf = 0.35	/*From paper*/
		 
svmat c1c2_results
keep if c1c2_results3!=.
keep c1c2_results*
drop if c1c2_results3==0
ren c1c2_results2 		award_rate_F
ren c1c2_results3		award_rate_M
ren c1c2_results1		c1_c2
sort c1_c2

gen delta_arm=abs(award_rate_M-$arm)
egen x=min(delta_arm)
su c1_c2 if delta_arm==x
gen c1_c2_norm=c1_c2-r(mean)

lab var c1_c2_norm "(Normalized) log(c{subscript:1}/c{subscript:2})"
lab var award_rate_M "Men's award CDF"
lab var award_rate_F "Women's award CDF"

line award_rate_M award_rate_F c1_c2_norm, clp(solid dash) lcolor(blue red) ///
	yline($arm , lcolor(blue) lp(solid)) yline($arf , lcolor(red) lp(dash)) ///
	yline($pdf , lcolor(red) lp(dot) lwidth(2pt)) yline($pdm, lcolor(blue) ///
	lp(dot) lwidth(2pt)) ///
	text(.16 4.2 "P{subscript:D}(M)",size(vsmall)) text(.24 4.2 "P{subscript:D}(W)",size(vsmall)) ///
	text(.375 3.6 "Women's award rate",size(small)) ///
	text(.49 3.7 "Men's award rate",size(small)) graphr(c(white)) xtitle("(Normalized) log(c{subscript:1}/c{subscript:2})",size(small))	ylabel(0(.2)1) xlabel(-4(1)4) legend(pos(6) col(2) region(lc(none)))
graph export $figures\figure3.pdf,replace


gen d_F=abs(award_rate_F-$arf)
gen d_M=abs(award_rate_M-$arm)
egen mdf=min(d_F)
egen mdm=min(d_M)
su c1_c2 if d_F==mdf
scalar p1=r(mean)
su c1_c2 if d_M==mdm
scalar p2=r(mean)
scalar alpha_check=p1-p2
***Footnote 44
di alpha_check

log close
